Time-dependent natural orbitals and occupation numbers 
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We report equations of motion for the occupation numbers of natural spin orbitals and show 
that adiabatic extensions of common functionals employed in ground-state reduced-density-matrix- 
functional theory have the shortcoming of leading always to occupation numbers which are indepen- 
dent of time. We illustrate the exact time-dependence of the natural spin orbitals and occupation 
numbers for the case of electron-ion scattering and for atoms in strong laser fields. In the latter 
case, we observe strong variations of the occupation numbers in time. 

PACS numbers: 31.15.Ew, 31.70.Hq, 34.80.-i, 32.80.Qk 
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The description of quantum many-body systems out 
of equilibrium has become an important research topic 
in nuclear, plasma and condensed-matter physics. The 
common interest of different fields in non-equilibrium 
quantum evolution is mainly driven by experimental and 
technological progress which raises questions such as how 
many-body systems evolve on transient or non-adiabatic 
time-scales, how they thermalize or which kind of trans- 
port phenomena are to be expected in extreme environ- 
ments. 

The Kadanoff-Baym equations [TJ|5J|3] provide a rigor- 
ous basis to investigate the dynamics of non-equilibrium 
many-body systems. Although the equations are known 
for some decades, so far only ab-initio solutions for very 
small atomic or molecular systems under special assump- 
tions have been achieved |4l |5] . For a more comprehen- 
sive ab-initio treatment of non-equilibrium systems with 
many degrees of freedom they are still out of scope of 
present day computing facilities. 

An alternative for the study of non-equilibrium pro- 
cesses in many-body systems is provided by time- 
dependent density-functional theory (TDDFT) ^. 
TDDFT is currently the method of choice for electronic 
systems out of equilibrium because it generally achieves 
decent accuracy at affordable computational cost. All 
physical observables are in principle functionals of the 
time-dependent density [71 [8]. However, in practice it 
is sometimes rather cumbersome to find approximate ex- 
pressions for observables of interest. Prominent examples 
are double and above-threshold ionization of atoms and 
molecules in strong laser fields, where no accurate func- 
tionals for the observables are known [9]. 

Recently, there has been renewed interest in reduced- 
density-matrix-functional theory (RDMFT) ^ [TT]- 
RDMFT is a promising candidate to treat long standing 
problems present in traditional density-functional theory 
(DFT). RDMFT describes very accurately the cohesion 
and dissociation of diatomic molecules which presents a 
difficult problem for DFT methods due to the importance 
of static correlation [T^] . The theory has also been suc- 



cessfully applied to the calculation of the fundamental 
gap [13]. 

In this work we consider the description of non- 
equilibrium many-particle systems in terms of time- 
dependent reduced-density matrices. These density ma- 
trices have the potential to capture the physics of strong 
correlations p4j in the time-dependent case and, further- 
more, observables can be more easily extracted from den- 
sity matrices than from the density alone. 

The time-evolution of the reduced density matrices 
is given by the Bogoliubov-Born-Green-Kirkwood-Yvon 
(BBGKY) hierarchy of equations of motion [TF. The 
first equation of the hierarchy has the form 



idaiiW; t) = (/io(l; t) ~ /io(l'; i)) 71(11'; t) 
+ f Kc(12) - «ec(l'2)) 72(121'2'; t)\2' 



:d2, 



(1) 



where the bare single particle Hamiltonian is given 
by ft,o(l;i) = — V^/2 -I- Wcxt(l) and Wcc(12) denotes the 
particle-particle interaction. The coordinates are written 
as combined space-spin variables 1 = (ri,(Ti) and we 
use J dl = ^^ J dJ^i- 1"^ general, similar equations in 
the BBGKY hierarchy connect the time-evolution of the 
N-th order reduced density matrix 7jv to the density 
matrix of order N-l-1. Since 2N cordinates appear in the 
equation of order N, the propagation of the complete 
hierarchy is even more involved than solving the un- 
derlying time-dependent Schrodinger equation (TDSE). 
Hence, truncations of the hierarchy are performed in 
practice. To close the system of evolution equations 
at a given order N it is therefore necessary to express 
the matrix of order N-l-1 in terms of matrices with 
order less or equal to N. This is the central point where 
approximations enter in a time-dependent reduced den- 
sity matrix theory. Note, that in the case of two-body 
interactions Wco(12) it is sufficient for the calculation 
of the ground-state energy to know the diagonal of the 
two-body reduced density matrix 72(1212;t — 0). As 
can be seen directly from Eq. ([l]), this is in contrast to 



the time-dependent case, where also ofF-diagonal matrix 
elements 72(121'2; t > 0) enter the equation of motion. 

By diagonalizing the one-body reduced-density matrix 
7i(ll';t) at each fixed point in time one obtains eigen- 
vectors ipj{l;t) and eigenvalues nj(t) which we term, in 
analogy to static RDMFT, time-dependent natural or- 
bitals and time-dependent occupation numbers, respec- 
tively [IS] . The spectral representation of the one-body 
matrix then takes the form 
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{n';t) = J2MtWk{'^;tM{i';t). 



(2) 



Since the eigenvectors (pj{l; t) form a complete set at each 
point in time, we can also express the two-body reduced 
matrix in the basis of natural orbitals of the one-body 
matrix 



72(12, l'2';i) = Y.j2.Mt)^ 

ijkl 



(3) 



ip,il;t)ip,i2;tM{l';t)ip^{2';t). 



With Eq. (I2]) the time-dependent natural orbitals and oc- 
cupation numbers have been introduced by a diagonal- 
ization procedure. It is now interesting to observe that 
the occupation numbers obey evolution equations which 
have the following form 

« ^k {t) = ^ l2,ijki (t) { ij I Wcc I kl ) (t) - c.c. (4) 

ijl 

where we have used ( ij \ Wee | kl ) (t) as shorthand for the 
coulomb integrals 



^,(1; <)^,(2; t)i;ee(12)^^(l; i)^r(2; t) dld2. (5) 



Next, we separate the coefficients ^2,ijki{t) into a mean- 
field and a cummulant part X2,ijki(t) |17j 

l2Ajki{t) = ni{t)nj{t){6ikSji ~ Sii6jk) + >^2,ijki{t)- (6) 

Inserting mj into the equation of motion for the occupa- 
tion numbers Q one finds, that Hartree and exchange- 
like contributions cancel. Consequently, we obtain that 
only the remaining cumulant part of the two-body re- 
duced density matrix determines the time-evolution of 
the occupation numbers 



i^kit) = ^X2A]ki{t){i3\vco\kl){t) 

ijl 



(J) 



For the closure of the BBGKY hierarchy up to level 
N=l an approximation of the two-body reduced density 
matrix in terms of the one-body matrix is required. One 
might be tempted to extend common functionals used 
in static RDMFT to the time-dependent case in a spirit 
similar to the adiabatic local density approximation of 



TDDFT [IS], i.e. for a slowly varying time-dependence 
in the Hamiltonian the ground-state functional is evalu- 
ated at the time-dependent density. In the case of time- 
dependent reduced density matrices this corresponds to 
replacing the ground-state occupation numbers by their 
time-dependent counterparts. Although this replacement 
is rather ad hoc, the hope would be that such functionals 
perform also reasonably well in non-adiabatic situations. 
The functional form of most commonly used ground- 
state functionals in RDMFT, written in the basis of the 
natural orbitals, can be summarized with the following 
expression 



l2,ijki — fijkiSikSji — g-ijkiSiiSjk, 



(8) 



which contains Hartree (dikSji) and exchange-like (SuSjk) 
terms. As example, for the Miiller functional [T^] we 
have ftjki = l/2ninj, gijki = l/2^n^nj, the self- 
interaction corrected functional of Goedecker and Um- 
rigar [TU] reads fijki = 1/2 {nirij - n^SijSikSa), giju = 



1/2 (V^ 



iT-iSijSikSii) and the BBCl functional of 



Baerends et al. [15^ has the form fijki — 1/2 n^ Uj , gijki — 
^n-njil/2 - Su6 -kil - %)e(z -N- e)e{j -N^ e)), 
where & denotes the usual Heaviside step function and 
< e < 1. In a similar fashion the BBC2/BBC3 func- 
tionals can be written in the form of Eq. (l8|. Note, that 
all functionals have the symmetry gijki — gjuk and all 
matrices fijki, gijki £^re real valued. 

Replacing the static occupation numbers which appear 
in Eq. dsl by their time-dependent counterparts and in- 
serting this approximation for the two-body matrix into 
the equation of motion for the occupation numbers Q 
we obtain 



ifikit) = Y.^fkm{t) - rt.jkimkj I I'ec I kj){t) 

+ ^(9]kk]{t) - 9jkkj{t)){jk I Wcc I kj){t), 



(9) 



which shows that all functionals of the form dsl) with real 
valued matrices fijki, gijki cause a zero right hand side 
in ([9]). Hence, if this class of approximations is used for 
the time-evolution of the one-body matrix 71 , the occu- 
pation numbers stay constant in time. This is a severe 
shortcoming of an adiabatic extension of present func- 
tionals of static RDMFT which needs to be addressed 
in the development of future functionals. Possible func- 
tional forms that lead to a non- vanishing right hand side 
in (|9| would be j2.ijki{t) = g{fik{t) fjiit) - fnit) fjkit)) 
or l2,i]ki{t) = giifijit) - fji{t)){fki{t) - fik{t))), where 
fij is a non-diagonal real-valued matrix and g some 
Taylor-expandable function. Alternatively, functionals 
with complex valued matrices could be employed. 

How do occupation numbers actually change in time? 
In the discussion above we have seen that present func- 
tionals of static RDMFT would lead to time-independent 
occupation numbers. Is this a sensible approximation? 
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FIG. 1: (a) Correlation entropy of e-He^ scattering for differ- 
ent interaction strengtlis A. (b), (c) Occupation numbers for 
the largest natural orbitals, respectively. 



To date the time-dependence of occupation numbers and 
natural orbitals is totally unknown. Take as an example 
the He atom and consider a transition from the ground 
state (where the two Is natural spin orbitals have occupa- 
tion numbers close to 1) to the first excited state (where 
the Is and the 2p spin orbitals have spin-degenerate oc- 
cupation numbers close to 1/2) by virtue of an optimized 
laser pulse. How do occupation numbers and natural 
orbitals evolve during the transition? At first sight two 
extreme scenarios are conceivable for this case: (i) the oc- 
cupation numbers stay nearly constant in time and only 
the orbitals change as follows: Is -^ Is, Is -^ 2p; or the 
opposite scenario: (ii) the orbitals stay close to the ini- 
tial orbitals and only the occupation numbers change to 
achieve the (ls,2p) final state. Particularly interesting 
is the fact that this is a transition from a weakly cor- 
related to a strongly correlated state (the latter being a 
superposition of two determinants) . 

To investigate the exact time-dependence of natural 
orbitals and occupation numbers we consider two pro- 
totypical cases: (i) e-He"*" scattering and (ii) the above 
mentioned transition of the Helium atom from the ground 
state to the first excited state induced by an optimized 
short laser pulse. In both cases we limit ourselves to 
reduced dimensionality so that the full time-dependent 
many-body Schrodinger equation can be solved numeri- 
cally. From the exact time-dependent many-body wave 
function we compute the one-body reduced density ma- 
trix and by diagonalizing the matrix at each point in time 
we obtain the natural orbitals and occupation numbers. 
We evolve both systems (i) and (ii) under a Soft-Coulomb 



(a) 75 

50- 
25 

S. " '~ 
-25 
-50 

-75 - 



(b) 


75 




50 




25 






ro 





X 






-2b 




-50 




-75 



50 100 150 200 250 300 

1 (a.yj 



50 100 150 200 250 300 
t[a.u.) 



FIG. 2: Space-time plot of the electron density for e - He"* 
scattering, (a) A = 1.0, (b) A = 2.0. 



Hamiltonian [2T] of the form 



m = 



pI , pI 2 



v/^f+I ^/x[+l 



(10) 



^(xi-X2)2 + r 

To vary the degree of correlation in the many-body 
wave function we introduce a coupling constant A in the 
Hamiltonian which controls the strength of the electron- 
electron interaction. 

Ref. [25j introduced the correlation entropy 



s{{nk}) :=-— rr7iln7i = -^ ^nfelnn^ 



(11) 



as a measure for the correlation which is contained in 
the one-body reduced density matrix. The correlation 
entropy s defined in this way is identical to zero for non- 
interacting particles and grows with increasing correla- 
tion in the system. In our time-propagations we calcu- 
late s({nfc(i)}) as function of time to monitor the degree 
of correlation which is present in the many-body wave 
function at a given instant in time. 

Case (i) - electron-ion scattering. For the initial state 
of the propagation we consider an antisymmetric spin sin- 
glet product wave function formed from a Gaussian wave 
packet ip{x) = exp( — (x — xq)'^ /a'^) exp{—kox) and the 
ground state (j)o{x) of the He+-ion. For all calculations 
we place the packet initially at a distance of xq = 15 a.u. 
away from the ionic core and give the scattering electron 
a momentum of kg = 0.3 a.u. which is pointing towards 
the ion. In Fig. [l] we plot the occupation numbers and 
the correlation entropy as function of time for different 
values of the electron-electron interaction strength A. Af- 
ter about i = 10 a.u. the wave packet has approached the 
ionic core and is passing the atomic nucleus. During this 
time the degree of correlation in the wave function is en- 
hanced. After the collision the transmitted and reflected 
waves are leaving the ionic center (cf. Fig. [2]) and the 
correlation entropy starts to decrease. As indicated by 
the correlation entropy for long times after the scattering 




FIG. 3: (a) Optimal laser pulse for the transition of the He- 
lium atom from the ground state to the first excited singlet 
state. The overlap of the propagated wave-function and the 
target state was 98.59% after 12 OCT iterations, (b) The two 
largest occupation numbers and (c) the correlation entropy as 
function of time. 



event, the many-body wave- function is again well repre- 
sented by a single Slater determinant. The occupation 
numbers deviate most strongly from their determinantal 
values (0 or 1) for A = 1.5. For A = 2.0 the electron- 
electron repulsion is already so strong that the incident 
wave packet is mainly reflected back. 

Case(ii) - optimal control. To study laser induced tran- 
sitions in the Helium atom we add an external dipole 
laser field of the form y(a;i,a;2) = (^i + X2) E{t) to the 
Hamiltonian H\=i. We use standard optimal control 
theory (OCT) [20] to find the optimal laser pulse with 
amplitude E(t) which drives the atom in a finite time- 
interval [0, iond] from the ground state to the first excited 
state. The full solution of the time-dependent many- 
body Schrodinger equation shows, that the actual transi- 
tion is a mixture of both anticipated scenarios: the nat- 
ural orbitals as well as the occupation numbers change 
as function of time during the transition. The occupa- 
tion numbers undergo large changes which refiects the 
multi-reference nature of the first excited state while the 
orbitals nicely reflect the quiver motion of the electrons 
in the laser field. This is displayed in Fig. |3] and Fig. |4] 
where we plot optimal laser pulse, correlation entropy, 
occupation numbers and the orbital density of the natu- 
ral orbital with the largest occupation number. 

In summary, we have presented equations of mo- 
tion for the occupation numbers of the natural spin or- 
bitals. We have shown that the adiabatic extension of 
present ground-state functionals of RDMFT to the time- 
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FIG. 4: Time-evolution of the orbital density of the first nat- 
ural orbital of Helium during the excitation from the ground 
state to the first excited state. The quiver motion of the elec- 
trons is nicely refiected in the density profile of the natural 
orbital. 



dependent domain yields always occupation numbers 
which stay constant in time. The exact time-evolution 
of the natural orbitals and occupation numbers has been 
illustrated for electron-ion scattering and for the Helium 
ls-2p transition. The exact analysis shows that sizable 
changes in the occupation numbers can occur during the 
time-evolution of the system. Approximations beyond 
the ones used for static RDMFT will therefore be necce- 
sary to reasonably capture the time-evolution of the one- 
body reduced density matrix. 
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